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Abstract 

An embedding method for solving the time-dependent Schrodinger equation is developed using 
the Dirac-Frenkel variational principle. Embedding allows the time-evolution of the wavefunction to 
be calculated explicitly in a limited region of space, the region of physical interest, the embedding 
potential ensuring that the wavefunction satisfies the correct boundary conditions for matching 
on to the rest of the system. This is applied to a study of the excitation of electrons at a metal 
surface, represented by a one-dimensional model potential for Cu(lll). Time-dependent embedding 
potentials are derived for replacing the bulk substrate, and the image potential and vacuum region 
outside the surface, so that the calculation of electron excitation by a surface perturbation can be 
restricted to the surface itself. The excitation of the Shockley surface state and a continuum bulk 
state is studied, and the time-structure of the resulting currents analysed. Non- linear effects and 
the time taken for the current to arrive outside the surface are discussed. The method shows a clear 
distinction between emission from the localized surface state, where the charge is steadily depleted, 
and the extended continuum state where the current emitted into the vacuum is compensated by 
current approaching the surface from the bulk. 

PACS numbers: 71.15.-m, 73.20.-r, 74.25. Jb, 78.47.J- 
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I. INTRODUCTION 



The question "How long does it take?" arises in many areas of condensed matter physics, 
in photoemission from surfaced, electron tunnelling through barrierd^, relaxation after 
the creation of a core hol^ and so on. In this paper we shall study the time-dependence 
of electron emission from a metal surface, using an embedding formalism which we have 
recently developecP and now further improved. 

The development of attosecond streaking spectroscopy, using ultra-short laser pulses, 
has raised interesting questions about the time it takes for electrons to be emitted in 
photoemissiorP, and recent computational studies have tackled this problemP^. In these 
time- dependent calculations, the surface and substrate are represented by a slab with a 
finite region of vacuum outside, and inevitably an electron excited by the time-dependent 
perturbation reaches the edge. How to deal with this is a familiar problem in time-dependent 
calculations in atomic physics, and it is often handled (as in the calculations just mentioned) 
by introducing an absorbing potential at the boundaries. It is important to eliminate spuri- 
ous back reflections, and exact boundary potentials have been developed which ensure that 
the wavefunction has the correct boundary conditions to match on to the rest of spacd^*^. 
This is in fact the idea of embeddin^^*^ we solve the Schrodinger equation in the region 
of space of interest, which we call region I, with an embedding potential at the boundary 
of the region to ensure that the wavefunction matches correctly on to the external region, 
region II. 

The original time-independent embedding methocP^ has proved useful for calculations 
of surface electronic structur^^ - only the surface region, the top few atomic layers plus 
the near-surface vacuum region are treated explicitly, the bulk crystal and the semi-infinite 
vacuum being replaced by embedding potentials derived from the substrate Green functions. 
The embedding method treats the bulk continuum states at the surface, as well as localized 
surface states correctly, and gives a very accurate description of surface electronic structur^^^. 
Our time-dependent embedding method builds on time-independent embedding and in the 
previous papei'^ we transformed the embedding potential as a function of energy into a 
function of time, and developed an embedded time-dependent Schrodinger equation for 
region I. 

This time-dependent embedding used the results of earlier papers in which the problem 
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of the exact termination of the spatially discretized Schrodinger equation was studiecP^li^. 
However, in our method, we use a basis set expansion of the wavefunction in region I, 
with time-varying coefficients, to solve the Schrodinger equation. Recent time-dependent 
calculations on electron transport through molecules have used self-energies to describe the 
coupling of the molecular wavefunctions to the metal electrode^. The self-energy is in fact 
an embedding potential in a tight-binding representatiorJ^, and this approach is related to 
our time-dependent embedding. 

We have made several advances in the time-dependent formalism since the earlier papei'^, 
and we present these here, together with results on the emission of electrons from bulk and 
surface states at the Cu(lll) surface. In section [ll| we shall give a new derivation of the 
method based on the Dirac-Frenkel variational principld^^n^, appropriate to solutions of the 
time-dependent Schrodinger equation. This avoids some of the weaknesses of our earlier 
derivation, which was based partly on analogy with time-independent embedding. In our 
application to emission of electrons from the Cu(lll) surface, we use a one-dimensional 
model potential due to Chulkov et al^^^ to describe the bulk, the surface, and the image 
potential outside the surface - this gives a good description of the electronic structure around 
the s-p band gap, with the Shockley and image potential-induced surface states, and has been 
widely used in surface calculation^^'^. In section III we shall calculate the time- dependent 



embedding potentials for replacing the semi-infinite bulk and vacuum regions, making use 
of some recent results for the time-independent embedding potentials appropriate to the 
Chulkov surface potentia P^*^° l Tests of the method, in which we apply a time-dependent 



surface perturbation, are described in section IV A crucial test is provided by the continuity 



equation - the number of electrons in the surface region plus the time-integrated currents 
across the boundaries of the embedded region should be constant. This is satisfied to a high 
degree of accuracy. 

In section |V] we give results for the electron current through the boundaries of the embed- 
ded surface region, into the bulk and into the vacuum, when a time- dependent perturbing 
potential is applied at the surface. We find that after initial transients, the current into 
the vacuum settles down, with an average value given by Fermi's Golden Ruld^. With a 
low frequency of surface perturbation, non-linear oscillations appear in the current, due to 
the interference between first order and second order processes. We also study the time it 
takes for the electrons to arrive at various distances from the surface, comparing the arrival 



time with the classical value, and for the time for the current to peak. Our most surprising 
results concern the comparison of currents into the bulk and into the vacuum. In the case 
of emission from a bulk continuum state, the current leaving the surface region into the 
vacuum is balanced by the current entering the surface from the bulk, so that the charge 
in the surface region stays constant on average; however, in emission from a surface state, 
currents of very similar magnitude leave the surface in both directions. Finally, in section 



VT] we give a brief outlook for the method. 

Atomic units are used throughout this paper, with e = h = irif. = 1. The atomic unit of 
time = 2.418884 x 10"^^ s, so that 1 fs = 41.34138 a.u. 



II. TIME-DEPENDENT EMBEDDING FORMALISM 




FIG. 1. Region I is embedded on to region II over surface S. The Schrodinger equation is solved 
explicitly in region I, region II being replaced by the embedding potential over S. 



The idea of embedding is that we solve the Schrodinger equation explicitly only in re- 
gion I, the rest of the system, which we call region II, being replaced by an embedding 
potential over the interface S (figure [T]). We developed the embedding method for the 
time-independent Schrodinger equation using a variational principle, and showed that the 
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wavefunction satisfies the following equation in region I (we include S within region I), 
1. 



2 V + V{r) ] 0(r) + 6{r - rg) 



2 9ns 



+ ^ dr's ^(rs, r's; e) + {E - e) ^) 0(r'^) 



E0(r). (1) 



Gq ^(rs, r'g; e) is the embedding potential, a generalized logarithmic derivative of the solution 
of the Schrodinger equation in region II at energy e, 



dri: 



s 



-2 / rfr'sGo^(rs,r's;e)^(r's), (2) 
Js 



where we take the normal derivative on S outwards from region I into II. The two terms 
in the round brackets in ([T|, the embedding potential evaluated at energy parameter e and 
the energy derivative term, g energy E, the energy at which we are solving the 

Schrodinger equation, to first order in [E — e). The surface terms inside the square brackets 
then vanish when matches in amplitude and derivative on to the solution of the Schrodinger 
equation in region II, giving the correctly embedded solution of the Schrodinger equation. 

We now turn to the time-dependent embedding problem, where we shall assume that the 
time- dependent potential, such as an external perturbation, is confined to region I. First 
we need a relationship analogous to ^ for the solution of the time-dependent Schrodinger 
equation in II. This is given 

where the time-dependent embedding potential is given by, 

I Go-^(rs,r's;6) 



Go^(rs,r's;t) = - / &exp(-iet) " ^ ^; ^ (4) 
^TT -le 

Equation ^ is almost the Fourier transform of (|2]), but it contains the time-derivative 
dijj/dt'; this compensates for the convergence factor of — 1/ie in the transform Q of the 
embedding potential. When we evaluate this transform, we take e to lie just above the real 
axis, above the singularities in the Green function, so that GQ^{t) = for t < 0. 

With the time-dependent version of the generalized logarithmic derivative (|3|, we can 
now derive the time-dependent embedded Schrodinger equation, and unlike our original 
derivatiorP, here we shall use a variational method. There are, in fact, several stationary 
principles for the time-dependent Schrodinger equation!^, and we choose the Dirac-Frenkel 



variational principl d^^*^^^ , which gives the functional variation 

d 



51= /rfr5^*(r,t) 



dt 



^(r,t) = 0, 



(5) 



for small variations of the time-dependent wavefunction. The Hamiltonian if in (jsj is 
given by 

H = -]^V^ + V{v,t), (6) 

and the integral is over the combined regions I + II. 

We use the same form of wavefunction \1/ as in the original embedding methocP^ - in 
region I, \1/ is given by a trial wavefunction (f){r,t), and in region II by the exact solution of 
the time-dependent Schrodinger equation ilj{r,t), which matches in amplitude on to 0(rs,t) 
over the boundary S between the two regions (0 and are assumed to be zero over the outer 
boundary of the combined region I + II, figure [T|. Evaluating the integrand in ([s]) with this 
form of the volume integral splits into the two regions, with an additional contribution 
on S coming from the kinetic energy operator — ^V^ acting on the discontinuity in derivative 
between and ip, 



51 



dv5<j)* 



H(j)-i 



.dcj)' 



dt 



+ / dr5i)* 
'ii 



Hi) -I 



.dip 



dt 



d(f) dip 
driQ drifi 



(7) 



But as ip is an exact solution of the time-dependent Schrodinger equation in region II, the 
second integral is zero. Also using ^ and the constraint that ip{rs,t) = (p{vs,t), the term 
involving dip/dris can be rewritten as 
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drsl 



dip 



s / "^s 

s 



f(rs,t)Go-(rs,4t-t')^^§^. (8) 



This assumes that on S is zero in the distant past, and ip itself is zero - in other words, 
ip grows into region II from region I as time proceeds. The variation of the functional then 
becomes an expression only involving (p in region I and on S, with the embedding potential 
accounting for region II, 



51= dr5(p*{r,t) 



H-i 



d_ 
dt 



[r,t) + 



'(rs,t) 



dns 



+ [drs [ dr's [ 5<P*irs,t)Go\rs,r's;t-t 



>,d< Pir's,t') 
dt' ' 



(9) 



We see that 51 = for arbitrary variations of 5(p* if (p satisfies the time-dependent 



Schrodinger equation in region I analogous to ([T]), 



V^ + \/(r,t) U(r,t) + 5(r-r5) 



1 d(j) 
2d^s 



Js J-oo 



dt' 



dt 



(10) 



We can also use (|9| to obtain the time-dependent matrix equation, with which we solve 
the embedded time-dependent Schrodinger equation in practice. Let us expand the trial 
function in terms of a set of basis functions Xii 



(11) 



- it is convenient both in the formalism and in coding if the basis functions are orthonormal 
when integrated over region I. Substituting this expansion into Q gives 

5/ = 5^KW^..(^KW-iE^<W^ + E^<W f dt%,{t-t')^, (12) 
where the Hamiltonian matrix Hij includes the surface derivative term in ([9 

Hijit) = ^ jdrVxiir) ■ Vxi(r) + J^dvXi{r)V{r,t)xjir) 
and the embedding matrix Sj^ is given by 



(13) 



Sij(t) = / drs / rfrsX»(rs)Go^(rs,rs;t)xj(rs). 
Js Js 



(14) 



Then SI = for arbitrary 6a* (t), and the variational principle is satisfied, when 



.dtti 



(15) 



- the time-dependent embedded Schrodinger equation in matrix form. This equation was 
derived in our earlier papei'^ in a less straightforward way. 

This formalism, which handles the evolution of a state initially localized in region I, needs 
modifying to study the time evolution of bulk continuum states at the surface, as well as 
the discrete surface states. As in earlier surface embedding worli^, we take region I to be 
the surface region, that is, the top few atomic layers plus the adjacent image potential -|- 
vacuum region, and region II to be the rest of the system, the semi-infinite crystal substrate 
on one side, and the semi-infinite vacuum on the other. The initial wavefunction then 



extends beyond region I - obviously so for continuum states, but even the exponentially 
decaying surface states extend some distance into the bulk and vacuum into region II. The 
time-evolution of extended states was also treated in our previous papei'^, but here we use 
the Dirac-Frenkel variational principle (jsj) to develop the formalism more consistently. 
The unperturbed wavefunction S(r,t) is a stationary state of the Hamiltonian Hq, 

Ho = -lv' + Vir), (16) 

where V{r) is the time-independent potential in I + II. We have 

E{r,t) = e(r) exp{-iEt), with Ho^ = (17) 

Let us now switch on the perturbing potential 6V{r, t) in region I at time t = 0, and see how 
this wavefunction develops. As before we shall write the full time-dependent Hamiltonian 
as H , and the evolving wavefunction as with 

H = Ho, ^(r,t)=H(r,t), t < 0, (18) 

and for t > 0, 



H = Hq + 5V{Y,t), ^(r,t) = S(r,t) + 0(r,t) r in region I 

H = Ho, ^(r,t) = E:(r,t) + V^(r,t) r in region II. (19) 

Again, 0(r, t) is the trial function in region I, and ip{r,t) is the exact solution of the time- 
dependent Schrodinger equation in region II, which matches on to 0(rs, t) over S. So S + ■?/' 
satisfies the time-dependent Schrodinger equation in region II, and the functional variation 
([5| simplifies to 



SI = / dr6^*{r,t) 



1 /■ , e,*. ,^<90(rs,^) 
+ o / drsS(j) (rs,t)- 



+ [drs /dr's /'',;0-(rs,i)Go'(rs,r's;i-i') ^'^3^;''' . OO (20) 

•/ S S •-'0 

- the same as ([9]), apart from the extra inhomogeneous term SVE in the first integral, and 
in the time integral we have taken a lower limit of t = 0, the time at which the perturbation 
is switched on. 
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We can finally write down the matrix equation for the time evolution of extended states. 



Substituting expansion (11) into (20) gives 



51 = Y,6a*{t)H,,{t)a,{t) + ^ 5a*(t)e,(t) - iJ^KW 

ij i i 

IJ 

where the Hamiltonian matrix now includes the perturbing potential, 



(21) 



H,,{t) = - jdvVxi{^) ■ Vx,(r) + jdvx,{v)[V{v) + 8V {v ,t)]xj{v) 
and vector Cj comes from the unperturbed state, 



(22) 



ei{t) = jdrx^{r)5V{v,t)'E{v,t). 



(23) 



The variational principle then gives the inhomogeneous time-dependent matrix equation. 



+ e,(t)=i^, t>0, 



(24) 



which we can solve for aj(t) with the initial condition that aj(t = 0) = 0. Knowing the Oj 
gives us 0, hence \l/ - the time-evolving wavefunction in region I, embedded on to the rest 
of the system. 



III. SURFACE EMBEDDING POTENTIALS 



We shall apply this formalism to study the excitation of electron states at a metal surface 
under the influence of a time-dependent perturbing potential. In the previous papei'^ we 
applied this formalism to a jellium surface, but here we shall consider a more realistic model 
of the surface, one which can support surface states as well as the bulk continuum. This is 
the one- dimensional model potential V{z) developed by Chulkov et alJ^^^, which has the 
form, 

Aicos{27Tz/a), z<0 
-Aio - A20 + A2 cos{l3z), < z < zi 

-A10 + A3exp[-a{z - zi)], Zi < z < z-,^ 

-Aio + (exp[-A(z - Zim)] - l)/4:{z - zi^), z > z-,^. 



V{z) 



(25) 
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region with z < 


a = 3.94 


Ai = 0.18889 






< z < zi 


(3 = 2.9416 


Aio = -0.43713 


A2 = 0.15905 


A20 = 0.40729 


Zl < Z < Zim 


Q = 0.63650 


Zl = 1.33499 


^3 = -0.51975 




Z > Z[Ya 


A = 1.27300 


Zim = 2.10562 







TABLE I. Parameters used in the Chulkov potential for Cu(lll) for different regions near the 
surface, a, Ai, Aiq, A2 and /3 are tabulated by Chulkov et al^, and the other parameters are 
derived; all parameters are given in atomic units. 

Here a is the interlayer spacing, parameter Ai reproduces the width of the bulk energy 
gap, AiQ is the inner potential (a difference with Chulkov et al.'^^ is that we take the average 
potential in the bulk as the energy zero, rather than the vacuum potential), and A2 and /3 are 
fitted to surface state energies. The remaining parameters are found by the requirement that 
the potential and its derivative are continuous across all the boundaries. The electrons are of 
course moving in three dimensions, but the motion parallel to the surface is free electron-like, 
and we shall only consider motion in the 2;-direction. Figure [2] shows the resulting surface 
potential for Cu(lll), the surface we shall be studying in this paper, calculated with the 
parameters given in table [1} 

We now derive the time- dependent embedding potentials for embedding the surface region 
on to the substrates on either side (figure [2]), G^^{t) for embedding on to the semi-infinite 
crystal substrate at Zc-, and G~^{t) for embedding on to the vacuum at z^. We start with 
the embedding potential for replacing the semi-infinite crystal, first calculating the energy- 
dependent embedding potential G~^{e) for subsequent Fourier transforming according to 
Q. As we have shown elsewher^^Sl^ ^j^e first step is to integrate the Schrodinger equation 
at energy e through a bulk unit cell in each direction to give us two independent solutions, 
01 and 02- Taking the unit cell to lie between Za and Zc (figure |2]), the boundary conditions 
for determining 0i and 02 are as follows, 

01 (-2a) = 1, <P'i{.Za) = 0, integrate from Za to Zc 

02 (-2c) = 1, 02 (-2c) = 0, integrate from Zc to Za- (26) 
The corresponding wave- vector k of the bulk band structure is then given by the remarkably 



10 



-15 -10 -5 5 10 15 

z(a.u.) 



FIG. 2. One-dimensional potential for modelling the Cu(lll) surface. At the surface region 
is embedded on to the semi-infinite crystal potential, and at it is embedded on to the vacuum 
Coulomb tail. The Schrodinger equation is integrated between Za and Zc to find G^^(e). 



simple expression, 

cos(fca) = '^^^"-^ + '^^^""^ (27) 

- a related result has been given by ButtP^, and earlier by Kohn!^. Using the property that 
it is a logarithmic derivative, it can be showrP^that the embedding potential for embedding 
on to the semi-infinite crystal to the left of Zc is given by 

' ^ ' 2[exp(-ifca)-0i(2,)]' ^ ' 

where W is the Wronskian. For causality, the Schrodinger equation is solved at e with 
a positive imaginary part (which may be infinitesimal), and k is chosen to have negative 
imaginary part, corresponding to a wave travelling or decaying into the crystal to the left. 
Figure |3] shows the embedding potential calculated in this way for the Cu substrate in the 
[111] direction, evaluated at Zc = —10 a.u., with an imaginary part in the energy of 2.5 x 10~^ 
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FIG. 3. Embedding potential for Cu(lll) evaluated at Zc = —10 £1.11.5 ^ 

function of Ree. Ime 
is kept fixed at 2.5 x 10^^ a.u. Curves a and 6, real and imaginary parts of G^^(e). 

a.u. The singularities are at the band edges. 

The next stage in finding G~^{i) is to take the Fourier transform of G^^(e), with the 
extra factor of — 1/ie Q. At large |e|, the embedding potential has free-electron behaviour, 
proportional to |e|^/^, and we can ensure rapid convergence of the integral by subtracting 
this off and adding on the corresponding time-dependent potential, 

Gl\t) = ^ I^Je ^-^^^[G-\e) - G]\e)] + G]\t), (29) 

where the free-electron embedding potentials are given b ^^ * ^" * ^^ * 



, J-e/2 e < - , t<0 



-iJe/2 e > ^ > 



We evaluate the integral in (29) by discretization between finite energy limits; we smear out 
the band edge singularities in the integrand by taking e with a small positive imaginary part 
(as in figure [3]), except in exp(— iet). 
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FIG. 4. Time-dependent embedding potential for Cu(lll) evaluated at Zc = —10 a.u. Curves a 
and b, real and imaginary parts of [G'^^(t) — Gj^{t)]; curves c and d, real and imaginary parts of 

Figure [4] shows the the time-dependent embedding potential for Cu(lll) G'^^(t), with 
the free-electron term subtracted; the free-electron embedding potential Gj^{t), with its 
inverse square root singularity as t — )■ 0, is also plotted. The embedding plane is taken at 
Zc = —10 a.u. (figure [2]). In the evaluation of the integral, the energy is discretized with 
Ae = 1.25 X 10~^ a.u. between energy limits of ±50 a.u., with an energy broadening of 
Ime = 2.5 x 10~^ a.u. The accuracy of this procedure is shown by the embedding potential 
satisfying causality very accurately, with G~^{t) almost exactly zero for t < 0. There are 
in fact some tiny Gibbs oscillations just below t = (figure |4]). The structure above t = 
is interesting, showing several periodicities coming from the structure in G^^(e) (figure [s]). 
The fundamental period is At ^ 15.2 a.u., which in Fourier transform corresponds to an 
energy of 0.41 a.u. - this is close to the energy difference of 0.42 a.u. between the bottom 
of the band and the top of the first band gap. 
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We now turn to the embedding potential to replace the vacuum region outside the crystal, 
in which asymptotically the electron feels the Coulomb tail of the image potential, 



in terms of the Chulkov potential parameters (|25j). We find the energy-dependent embedding 
potential from the logarithmic derivative expression (g, which in one-dimension simplifies 

to 



where ip{z) is the outgoing or decaying solution of the Schrodinger equation. This is a com- 
bination of the regular and irregular Coulomb functions Fo and Go with angular momentum 
L = 0, in the notation of Abramowitz and StegurP^, 

^(z) = H^{r^,p) = Go{v.p)-iFo{v,p) (33) 

with arguments given by 

p = A/2(e + Aio)(2im - z) 

V = — , ^ (34) 
Thompson and Barnett^ give a rapidly converging continued fraction expression for 



Hq' / Hq , from which we can immediately obtain G~^{e) using (32). This is shown in 
figure |5| evaluated at = 10 a.u. with an image plane at = 2.1056 a.u. In this figure 
we measure energy relative to the vacuum zero, and the structure in G^^ just below zero 
comes from bound states of the Coulomb potential. 

To find the corresponding time-dependent embedding potential G~^(t), we use the same 



procedure as in (29), subtracting the free-electron embedding potential inside the Fourier 
transform integral, and adding the time-dependent free-electron embedding potential out- 
side. The Fourier transform in the case of the image potential is more demanding computa- 
tionally than for the crystal potential, and we use an energy discretization of Ae = 1 x 10~^ 
a.u. between energy limits of ±50 a.u., with an energy broadening of Ime = 2.5 x 10^"^ a.u. 
However, even with this tiny value of Ae, G~^{t) turns out to be finite, though small and 
constant, for t < 0. This error, which is proportional to Ae, is actually constant over the 
whole range of t, positive as well as negative. This means that we can improve the accuracy 
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FIG. 5. Vacuum embedding potential evaluated at Zy = 10 a.u., with = 2.1056 a.u., as a 
function of Ree measured relative to the vacuum zero. Ime is kept fixed at 2.5 x 10~^ a.u. Curves 
a and b, real and imaginary parts of G^^(e). 

- essentially going to the limit of Ae = - simply by subtracting off the constant value of 
Gy at negative t. The resulting time-dependent embedding potential is shown in figure [g} 
After the correction, the embedding potential is almost exactly zero for t < 0, with even 
smaller Gibbs oscillations than in figure |4} 



IV. TIME-DEPENDENT TESTS 

In this section we shall test the accuracy of the time-dependent embedding formalism, 
applying it to electron emission from the Cu(lll) surface. But first we solve the time- 
independent Schrodinger equation to find the static electronic structure of the surface, using 
the Chulkov potential and the energy-dependent embedding potentials G^^{e) and G~^{e) 
calculated in the last section. We determine the Green function G{z, z'\ e) satisfying the 
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FIG. 6. Time-dependent vacuum embedding potential evaluated at Zy = 10 a.u., with = 2.1056 
a.u. The image potential is measured from the vacuum zero. Curves a and b, real and imaginary 
parts of [G-Ht)-G^^{t)]. 

embedded Schrodinger equation at complex energy e with a small imaginary part, e = E+irj, 



2 dz 



2 +{V~e)G{z,z';e) + 6{z-z,) 



IdG 
'2lh 



+ G-\e)G{z,,z';e) 



+6{z - z^) 



IdG 
2lh 



+ G-\e)G{z,,z'-e) 



5{z — z'\ z, z' in region I. (35) 



The local density of states - the electron density at energy ii^ - is then given 

a{z,E) = —lmG{z,z;e). 
G{z, z'\ e) is expanded in terms of basis functions XiiAi 



(36) 



G{z,z'\e) = ^Gij{e)xiiz)xjiz') 



and then (35) becomes the matrix equation. 



[-f^ij + Sjj(e) — eSij] Gjk — 5ik, 



(37) 



(38) 
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with the Hamiltonian matrix given by 



Hj^ = — 
u 2 



dxi dxj 



/ dzXiiz)V{z)xjiz), 



(39) 



the embedding matrix by 



Sij(e) = G'^ {e)xiizc)Xjizc) + ie)xiizv)Xjizv 



and the overlap matrix by 



Sij 



We use trigonometric basis functions, 

Xm{z) -- 



dzxi{z)xj{z)- 



COS , m even 
sin m odd 



(40) 



(41) 



(42) 



where is measured from the mid-point of region I, 

>. H" 
C = Z : 



(43) 



and a value oi D > {z^ — z^) /2 to give a range of logarithmic derivatives at Zc and z^ for 
matching on to the embedding potentials. In the time-dependent applications we construct 



an orthonormal basis from (42), with the eigenvectors of the overlap matrix as coefficients. 



Integrating the local density of states (36) over the embedded region, we obtain the 



surface density of states for Cu(lll) shown in figure [7j In this calculation we take region I 
between Zc = —10 a.u. and z^ = +10 a.u., and use 40 basis functions, defined with D = 12 
a.u. in (42); the imaginary part of the energy, which broadens the discrete states, is 1 x 10~^ 



a.u. These results, which correspond to zero wave-vector parallel to the surface in three 
dimensions, show the Shockley surface state at E = 0.2415 a.u., close to the bottom of the 
band gap at 0.2201 a.u, with the first image state aX E = 0.4072 a.u., immediately below 
the top of the gap at 0.4087 a.u.l^ the higher members of the image series overlap the 
continuum and become resonances. For comparison, the Fermi energy is at 0.2556 a.u.l^ 
and the vacuum zero at 0.4371 a.u. (table [l]). The surface state energies are almost identical 
to the values given by Chulkov et alW^, and close to experimental result^^^^. A caveat is 
that although the Chulkov potential reproduces the electronic states near the s-p band gap, 
in particular the Shockley and image surface states, it does not describe the d-states, which 
lie 1 - 3 eV below the band gap^^, nor does it give the band width accurately. 
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FIG. 7. Surface density of states of Cu(lll), the local density of states integrated through region 
I, taken between Zc = —10 a.u. and = +10 a.u. Energy broadening is 1 x 10^^ a.u. 



We now return to solving the time-dependent embedded Schrodinger equation (24) for 
Cu(lll), with a surface perturbation given by 

SV{z,t) = Aexp-{z'^/E)sm{ut). (44) 

This perturbation is switched on at t = 0, and we follow the subsequent time-evolution 



of 4>{z^t) (19), with expansion coefficients a^, given by equation (24). The unperturbed 
wavefunction ^{z) ( 17p8 ), which enters the time-evolution equation through vector Cj (23), 
can be efficiently calculated using the Numerov methocP^, integrating inwards from the 



vacuum. The matrix elements of the Hamiltonian Hij(t) (22) are given by (39), with 6V{z, t) 
included, 



dz 



dXi dXj 
dz dz 



dzx^{z)[Viz) + SViz,t)]xJiz), 



and the embedding matrix (14) has the same form as (40) 



= {t)xt{zc)Xj{zc) + G„ {t)xi{zy)xj{zv) 



(45) 



(46) 
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To develop a simple time evolution algorithm, we write equation (24) in abstract notation, 



— = -i[Ha + T], where r = e+ / dt'i:{t - t') — . 
dt J Q dz 

Then to first order in the time interval 5t,and using central differences, we have 



a{t + 5t) = a{t) - iStH{t + 5t/2)a{t) - iStT{t + St/2). 



But to first order we can also write 



a{t + 5t) = a{t) - i5tH{t + 5t/2)a{t + 5t) - i5tT{t + 5t/2), 



(47) 



(48) 



(49) 



and adding (48) and (49) gives 
ait + 6t)= 1 + 1— Hi t + — 



1 



.5t-f 5t 
i—Hl t + — 
2 V 2 



a{t) -i5tTit + 



St 



(50) 



This is the formula we use to advance the wavefunction expansion coefficients forward in 
time. 



Although (50) turns out to be very accurate - much more so than the time evolution 



method we described in our earlier papeil^ - it is only stable for small 6t. Because the 



trigonometric functions in (42) are defined with respect to 2D, but only used within the 



smaller embedding range, there is the risk of overcompleteness, and the closer we are to 
linear dependency in the basis set, the smaller the time interval 6t needed for stability. 
Taking an embedding range {zc — z^) of 20 a.u., and 2D = 24 a.u., 6t = 0.002 a.u. gives 
stability with up to 40 basis functions. With an embedding range of 40 a.u. and 2D = 44 
a.u. this time interval gives stability with up to 75 basis functions. Bringing 2D closer to 
the embedding range helps matters, and we find that with {zc — z^) = 20 a.u., 2D = 22 a.u., 
and 30 basis functions we can go up to 6t = 0.016 a.u. In fact this is very accurate: figure 



8] shows the electron density |\l/(2;,t)p (19) of the Shockley surface state at E = 0.2415 a.u. 
(figure [T]), calculated with this latter basis set and the large time interval, at t = 200 a.u. 



after the application of the surface perturbation (44) with A = 0.2 a.u., S = 2 a.u. and 



frequency u = 0.5 a.u. We see that it is almost identical to the calculation with the same 
embedding range, 2D = 24 a.u., 40 basis functions, and 6t = 0.002 a.u.; the difference is 
only visible at the enlarged scale of the inset. It is tempting to take 2D = {zc — z^) to 
eliminate the possibility of overcompleteness. However, for a given size of basis set this is 
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FIG. 8. Shockley surface state electron density at t = 200 a.u. after applying the surface pertur- 
bation, calculated between Zc = —10 a.u. and z^j = +10 a.u.: solid line, 30 basis functions with 
2D = 22 a.u. and 6t = 0.016 a.u.; short-dashed line, 40 basis functions with 2D = 24 a.u. and 
6t = 0.002 a.u. Long-dashed line shows the unperturbed density. 

less accurate, and it is much better to use a basis set with 2D > (z^ — z^), giving a range of 
derivatives at the embedding surfaces. 

To test the embedding method itself, we compare the time evolution of states calculated 
with different embedding surfaces. First we consider the system with Zc = —10 a.u., z^ = 
+10 a.u., taking the basis set with 2D = 24 a.u. and 40 basis functions, and a time interval 
of 6t = 0.002 a.u. In the second case region I extends deeper into the crystal, with Zc = —20 
a.u., and further into the vacuum, with z^ = +20 a.u.; we use a basis set with 2D = 44 
a.u., 70 basis functions, and the same time interval of 6t = 0.002 a.u. Let us consider the 
time evolution of the continuum state with E = 0.1 a.u. (figure [?]), with the same surface 
perturbation as before, A = 0.2 a.u., S = 2 a.u., and u = 0.5 a.u. The results for the two 
calculations are shown in figure |8| and once again the difference between the two curves is 
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FIG. 9. Electron density of continuum state with E = 0.1 a.u. at t = 200 a.u. after applying the 
surface perturbation: solid line, = — 10 +10 a.u., 40 basis functions with 2D = 24 

a.u.; short-dashed line, Zc = —20 a.u., z^, = +20 a.u., 70 basis functions with 2D = 44 a.u. In 
both cases time evolution proceeds with St = 0.002 a.u. Long-dashed line shows the unperturbed 
density. 

only visible in the magnified inset. 

Another interesting test is to look at the number of electrons in the embedded region, 
Q{t), as a function of time, 

Q(t)= rdz\^!{z,t)\\ (51) 

J Zc 

Continuity of charge requires 

Q{t) + Jc{t) + Jv{t) = constant, (52) 
where Jc and J„ are the time-integrated currents leaving region I across the left- and right- 



21 



O 0.8 




FIG. 10. Solid line, Q{t), Shockley surface state electron number in the embedded region as a 
function of time; dashed line, [Q{t) + Jc(t) + Jvit)], electron number + time-integrated current 
across the embedding surfaces. Zc = —20 a.u., = +20 a.u., 70 basis functions with 2D = 44 a.u.; 
time evolution proceeds with 6t = 0.002 a.u. 



hand embedding surfaces, 



.In 



=Flm — 

oz 



(53) 



We use the embedding formula ^ to calculate the normal derivative in the expression for the 
current. Taking region I between ±20 a.u., 70 basis functions with 2D = 44 a.u., and using 
6t = 0.002 a.u., we calculate the time evolution of the surface state at 0.2415 a.u. under the 
influence of the same surface perturbation as in the previous tests. The results are shown in 



figure 10 , and we see the steady decrease in the surface state charge in the surface region as 



time progresses. However, the continuity equation (52) is satisfied remarkably accurately 



with [Q + Jc + Jv] staying almost perfectly constant - the maximum variation is less than 
10~^ across the whole time interval. 
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V. ELECTRON EMISSION RESULTS 



Having tested the embedding method, we turn to the physics of electron emission from 
the Cu(lll) surface. As in the previous section, we apply the surface perturbation given by 



(44) - this might correspond to the perturbation in surface photoemission with p-polarized 



light, when the E-field of the incident light has a large component normal to the surfac^^. 
Parameter S is taken as 2 a.u., and we shall use amplitudes A of 0.01 and 0.1 a.u. (for 
comparison, the amplitude of the pseudopotential in the Chulkov potential for Cu(lll) is 
0.19 a.u., and the surface barrier is 0.44 a.u. (figure [2])). The results are calculated using an 
embedded region between Zc = —20 a.u. and = +20 a.u., with 70 basis functions defined 
with 2D = 44 a.u. 
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FIG. 11. Emission from the continuum state at E = 0.1 a.u., showing current measured at = 20 
a.u. as a function of time. Curves a, u = 0.4 a.u.: soHd line, x 100 with A = 0.01 a.u.; dashed 
Une, with A = 0.1 a.u. Curve b, x 100 with uj = 0.8 a.u., A = 0.01 a.u. 



We begin with emission from the continuum state at = 0.1 a.u., normalized so that the 
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electron density equals the local density of states (36), that is, energy normalization. Figure 



11 shows the current j„ emitted into the vacuum across the right-hand embedding plane 
at at two frequencies of perturbation, u = 0.4 a.u. and 0.8 a.u.; we use amplitudes of 
A = 0.01 a.u. and 0.1 a.u. at w = 0.4 a.u., and 0.01 a.u. aX u = 0.8 a.u. First the results for 
small amplitude, shown by solid lines in the figure: there is a time-delay before the current 
reaches - we shall discuss this later - but after a large initial peak, the current shows 
oscillatory transients about a constant current. By curve-fitting (we fit a straight line to the 
time-integrated current J^), we estimate that the average current after the initial peak is 
1.0 X 10~^ a.u. for u = 0.4 a.u. and 2.65 x 10^^ a.u. for u = 0.8 a.u. We can also calculate 
the long-term average current for small perturbations from Fermi's Golden Ruld^^'^, 

J. = ^\{¥V\f)\^ (54) 

where \i) is the initial state and |/) the final state, the well-known (though not particularly 
obvious) time- reversed LEED state, kj is the wave- vector of the final state (the incident 
free-electron wave in the LEED state is exp{—ikfz)), coming from the density of final states 
in one dimension, l/{2'Kkf) - the 27r cancels with the prefactor of 27r in the Golden Rul^22_ 
The Golden Rule expression then gives J„ = 9.62 x 10~^ a.u. for u = 0.4 a.u., and jy = 
2.65x 10~^ a.u. for u = 0.8 a.u., both ai A = 0.01 a.u., in agreement with the time- dependent 
calculations. The decrease in current with increasing frequency is due to a combination of 



density of states and matrix element effects, as we can see from (54). 

The lower the excitation frequency, the longer the period of the transients. Increasing the 
perturbation amplitude by a factor of 10 at a; = 0.4 a.u., the current shows very pronounced 
oscillations of much shorter period on top of the long period oscillations (dashed line, figure 



11). In fact even with A = 0.01 a.u., these short-period oscillations are visible. Comparing 



the current for A = 0.01 a.u. multiplied by a factor of 100 with that for A = OA a.u. (curves 



a, figure 11), we see that the amplitude of the short-period oscillations is non-linear in the 
intensity of the perturbation, varying as the cube of the perturbation amplitude A. This 
shows that these oscillations are the result of interference between first-order and second- 
order excitations from the bulk state; the frequency of the short-period oscillations is exactly 
0.4 a.u., as we would expect for an interference term. This higher order term in the current 
is much less pronounced in emission with u = 0.8 a.u. 



From the curves for jv{t) shown in figure 11 it is clear that there is no unique time at 
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which the electrons start to arrive at Zy, where we measure the current - do we mean the 
time at which the current differs from zero, or perhaps when the first peak arrives? The 
first is ill-defined, as the increase from zero is gradual, and numerical limitations mean that 
we cannot determine it accurately. The second is too late, especially for the emission curves 
at cj = 0.4 a.u., with a final state energy just above the vacuum zero. However, we can 
define a convenient measure from the straight-line fit to Jv{t)-, which we used to find the 



average current. Figure 12 shows the time-integrated current Jvif) evaluated at the right- 
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FIG. 12. Time-integrated current at = 20 a.u. emitted from the continuum state &t E = 0.1 
a.u., uj = 0.8 a.u., A = 0.01 a.u. Solid line, Jvit); dashed line, linear fit beyond t = 80 a.u. 

hand embedding plane in emission from the continuum state with E = 0.1 a.u., excited 
by the surface perturbation with frequency u = 0.8 a.u. and amplitude A = 0.01 a.u., 
fitted by a straight line beyond t = 80 a.u. We define the effective arrival time ty for the 
current as the intercept of the linear fit with the time axis - in this case, = 18.3 a.u. This 
compares with the classical arrival time, evaluated from the velocity of the electrons in the 



final state, of 20.8 a.u. Figure 13 shows jy{t) measured at different values of Zy, together 

25 



with the arrival times determined from the hnear fits: we see that our effective arrival 
times lie close to the start of the steep initial rise in current. The classical arrival times 
are about 2 a.u. later, still not long after the start of the current rise, and well before the 




t (a.u.) 



FIG. 13. Current emitted from the contimmm state at E = 0.1 a.u. measured at different planes 
as a function of time, m = 0.8 a.u., A = 0.01 a.u. Curve a, z„ = 20 a.u.; b, = 30 a.u.; c, = 50 
a.u. The vertical lines show the effective arrival times ty from fits to the time-integrated current. 

first peak. The time-structure of jv{t), with its transients, presumably reflects the way that 



the perturbation (44) is switched on suddenly at t = 0: a spectrum of energies is excited 



until the system settles down. This is why the initial rise in current widens when jv{t) is 



measured at larger Zy (figure 13) - j„ disperses. 

Our most interesting results concern the electron number Q{t) in the embedded region, 
and how this is related to the time-integrated currents across the bulk and vacuum embed- 



ding surfaces - we have already studied [Q{t) + Jdt) + Jvif)] in section IV Figure 14 shows 
Q(t), Jc(t), and Jv(t) for excitation from the continuum state with energy = 0.1 a.u. by 
the surface perturbation with u = 0.8 a.u. and amplitude A = 0.1 a.u. (we use a larger 
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FIG. 14. Upper figure, electron number Q{t) in the embedded region in emission from the contin- 
uum state at E = 0.1 a.u., with uj = 0.8 a.u., A = 0.1 a.u. Lower figure, time-integrated currents 
across the embedding surfaces: solid Hne, Jv{t); dashed line, Jc{t). 



amplitude to emphasize the results). Both Q{t) and Jc{t) show continuing short-period os- 
cillations, but if we average these by eye we see that after a decaying transient, Q settles 
down to a fairly constant value. This is the result of the current leaving the surface region 
into the vacuum across being balanced by the current approaching the surface from the 
bulk across Zc - the slopes of Jv(t) and Jc(t) are almost equal and opposite (remember, 
positive current means charge leaving the surface region). So in emission from a continuum 
state, the wavefunction is replenished from the bulk, and the electron number at the surface 
reaches a steady state. 

The behaviour is quite different in emission from the Shockley surface state at E = 0.2415 



a.u. Figure 15 shows the number of electrons in the surface region and time-integrated 



currents, with uj = 0.6585 a.u. (giving the same final state energy as in figure 14) and 



A = 0.1 a.u. We see that Q drops steadily, as the surface state is depopulated by currents 
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FIG. 15. Upper figure, electron number Q{t) in the embedded region in emission from the surface 
state at E = 0.2415 a.u., witli lo = 0.6585 a.u., A = 0.1 a.u. Lower figure, time-integrated currents 
across the embedding surfaces: solid line, Jv{t); dashed line, Jc{t). 

leaving the surface into the bulk as well as into the vacuum: Jc{t) is roughly equal to Jv{t), 
both directed away from the surface. The currents correspond to electron momentum, so 
the near-equality of Jc{t) and J„(t) approximately conserves momentum (the surface step 
also provides a source of momentum'^^). On the other hand, in emission from the continuum 
state momentum is directed from the bulk into the vacuum. 

One of the pleasing features of our work is that the embedding method describes the 
different time-evolution of the surface state and the continuum state, within the one-electron 
picture we are using. But we must remember that in real life, a state is depopulated in one 
go by the emission of a single electron, and many-body effects repopulate the resulting 
hole. The finite mean-free path of the electrons also means that the bulk reservoir will not 
repopulate the continuum state as efficiently as our calculation suggests. The difference 
between a localized surface state and an extended continuum state may be more apparent 
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than real after all. 



VI. OUTLOOK 

Having demonstrated and applied the method, we are aware of its limitations. Apart 
from the fact that it neglects many-body effects, there are restrictions to the embedding 
method itself. Chief of these is the fact that the time for evolving the wavefunction through 
5t is proportional to t, because the embedding term in the Hamiltonian involves an integral 



with the upper limit t (24). It is also quite time-consuming to evaluate the time- dependent 
embedding potential for the image potential region, though we can certainly go well beyond 
the present time-limit of 200 a.u. (about 5 fs), even with a desk-top computer. To reach 
much longer times further developments of the method are necessary. In the meantime 
there is much to do with the existing surface programs, considering excitation by ultra-short 
light pulses, multi-photon processes and the like. Our time-dependent embedding method 
should have other applications, for example to adsorbates^ and metallic contact^, where 
time-independent embedding has proved useful. 
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